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Abstract. We elaborate on a general method that we recently introduced for characterizing the "natural" 
structures in complex physical systems via a multi-scale network based approach for the data mining 
of such structures. The approach is based on "community detection" wherein interacting particles are 
partitioned into "an ideal gas" of optimally decoupled groups of particles. Specifically, we construct a 
set of network representations ("replicas") of the physical system based on interatomic potentials and 
apply a multiscale clustering ("multiresolution community detection") analysis using information-based 
correlations among the replicas. Replicas may be (i) different representations of an identical static system 
or (ii) embody dynamics by when considering replicas to be time separated snapshots of the system (with a 
tunable time separation) or (iii) encode general correlations when different replicas correspond to different 
representations of the entire history of the system as it evolves in space-time. Inputs for our method are 
the inter-particle potentials or experimentally measured two (or higher order) particle correlations. We 
apply our method to computer simulations of a binary Kob- Andersen Lennard- Jones system in a mixture 
ratio of A80B20, a ternary model system with components "A", "B", and "C" in ratios of AssBrCs (as in 
AlssYrFes), and to atomic coordinates in a Zr8oPt2o system as gleaned by reverse Monte Carlo analysis of 
experimentally determined structure factors. We identify the dominant structures (disjoint or overlapping) 
and general length scales by analyzing extrema of the information theory measures. We speculate on 
possible links between (i) physical transitions or crossovers and (ii) changes in structures found by this 
method as well as phase transitions associated with the computational complexity of the community 
detection problem. We briefly also consider continuum approaches and discuss the shear penetration depth 
in elastic media; this length scale increases as the system becomes increasingly rigid. 

PACS. 89.75.Fb - 64.60.Cn - 89.65.-s 



1 Introduction 

This article constitutes a longer companion work to an 
earlier summary [1 in which the basic notions to be de- 
tailed here for complex physical systems were succinctly 
outlined. We begin by briefly reviewing a special class of 
complex physical systems that is of great fundamental 
and technological importance - that of amorphous ma- 
terials. From a practical standpoint, these materials often 
have industrial processing and preparation advantages 
|3] relative to crystalline systems enabling, e.g., greater 
solubility of pharmaceuticals [4] and many other advan- 
tages ^|5j. Below, we list several specific complex amor- 
phous systems, (i) Metallic glasses can be stronger than 
their respective crystalline structure and exhibit interest- 
ing electrical, chemical, and magnetic properties [3]. (ii) 

^ Corresponding author: zohar@wuphys.wustl.edu 



Phosphate glasses are of great use in biomedical applica- 
tions and chalcogenide glasses are of vital importance in 
optical recording media such as Blu Ray technology [6 . 
(iii) Far more recent and exotic challenges involve incom- 
mensurate complex electronic structures found in systems 
such as the high temperature superconductors [7^8 . The 
understanding of the character of such non-trivial struc- 
tures is a problem of considerable interest in disparate 
arenas. 

In perfect crystals, the natural system scales are evi- 
dent by the regular ordering of the lattice. The fundamen- 
tal unit cells of a crystal typically involve several atoms 
that are replicated in a simple pattern to span the entire 
system. There are no intermediate scale structures within 
the system from the atomic scale of the lattice up to the 
complete single crystal. Identifying the basic periodic unit 
cells is vital to the understanding of all crystalline solids. 
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This simplicity of structure enables an understanding of 
crystalline solids in great detail. Early on, the existence 
of specific unit cell structures were postulated to exist in 
crystals based on the sharp facets and other macroscopic 
properties of large crystals. 

There are more complex systems in which new struc- 
tures appear on additional intermediate scales between the 
atomic-scale and the macro-scale of the system. In recent 
years, scientific exploration has endeavored to understand 
a vast array of such complex materials that do not have a 
simple theoretical starting point. 

As alluded to above, some of the best known complex 
materials are glasses. Liquids that are rapidly cooled ("su- 
percooled") below their melting temperature avoid crys- 
tallization and instead become quenched into an amor- 
phous state. On supercooling, liquids may veer towards 
local low energy structures (that cannot, on their own, 
be globally replicated to fill space without the inclusion 
of other structures) before being quenched into an amor- 
phous state. Local low energy structures such as those 
formed by icosahedral packing are indeed observed in metal- 
lic glasses OHO . Due to the lack of a simple crystalline 
reference, the structure of glasses is extremely hard to 
quantify beyond local scales. 

Many theories of glasses rely on the hypothesis of nat- 
ural structures in the glass pT| [T2 l [T3 l [T4] . Actually finding 
such structures in a general way has been more elusive. 
How then may we detect and best characterize the most 
notable structures in amorphous systems? In the current 
work, we introduce a general framework to answer this 
question with specific applications to model glass form- 
ers. 

The outline of this article is as follows. Initially, we 
present some background information concerning the pur- 
suit of characterizing structures in glasses and the ba- 
sic features of our method in Sect. |2j Details of the pri- 
mary simulated systems are presented in Sect. [Sj The 
community detection and multi-scale community detec- 
tions methods are explained in Sect. [4j Our multireso- 
lution method is applied to physical models in Sect. [5] 
We summarize our findings in Sect. [6) Details regarding 
the applied information measures are given in Appendix 
|A]and those concerning our overlapping dynamics proce- 
dure in Appendix [B] we used this overlapping commu- 
nity detection to augment the community partitions de- 
termined by the methods in Sect. [4j Appendix [C] high- 
lights the simple (yet often overlooked) fact that prepeaks 
in the structure factor do not constitute a necessary con- 
dition for medium range order. In Appendix [Dj we dis- 
cuss how we handle "overlapping nodes"- nodes that are 
common to one or more "communities" within the com- 
munity detection methods that we employ. The remaining 
appendices elaborate on additional test cases and various 
facets of our method. We conclude by discussing struc- 
tures in space-time and their general properties. We dis- 
cuss the detection of multi-scale structure in space-time 
in Appendix [R| In Appendix [Sj we consider continuum 
approaches and discuss the divergence of a general length 



scale- the "shear penetration depth" - as a supercooled 
liquid becomes quenched into a rigid glass. 

2 Background 

2.1 Lightning summary of numerous current 
approaches 

Existing work in the pursuit of understanding the glass 
transition is vast, spanning decades and affecting many 
fields of science and engineering. Glass formers exhibit 
several unique common features [15 . Glasses demonstrate 
short range order (SRO) and medium range order (MRO) 
structures, but no easily discernible static long range or- 
der exists. A striking property of glass formers (especially 
of the so-called "fragile" glass formers) is that their re- 
laxation times (as measured, by, e.g., their viscosity) can 
increase by many orders of magnitude over a relatively 
narrow temperature range. This dramatic slowing down 
is not associated with the usual distinguishing measures 
of conventional thermodynamic phase transitions. These 
systems have rich energy landscapes with an exceptionally 
high number of metastable states [16 ,17,18 ,19,20,21,22 . 

Given the broad appearance of glass-related states, dif- 
ferent frameworks have been explored to work towards a 
"universal" characterization of the glass transition. Many 
theoretical approaches, e.g., p^, 23, 24, 25, 26, 27, 14, 13 have 
been developed over the years. The notable theory of ran- 
dom first order transitions (RFOT) investigates mosaics of 
local configurations [15",'23". RFOT is related [M] to theo- 
ries of locally preferred structures [14, 13,28, TT|[T2]. Other 
theories seek a similar measures of structure. Amongst 
many others, these approaches include spin glass type 
analysis [24], theories of topological defects and kinetic 
constraints [14, 13,29,30,31 , and numerous ingenious ap- 
proaches summarized in excellent reviews, e.g., [19,32j 
[33] . Formally, as demonstrated in [34 , a growing static 
length scale is associated with the diverging relaxation 
times in supercooled liquids. Some works indeed found 
indications of increasing correlation lengths (static and 
those describing dynamic inhomogeneities) as the temper- 
ature was lowered [35]l36] [37 | [38] . Static correlation lengths 
were amongst other approaches, notably, examined in terms 
of (i) "point-to-set" correlations [39[ l4Q] as well as (ii) pat- 
tern repetition [41 . We will briefly discuss these measures 
later on. 

In metallic glasses, early work to ascertain local struc- 
tural used a dense random packing model [42 . It was later 
established that such structures are better represented 
by an efficient cluster packing (ECP) model jj3l[44l[45] . 
SRO features were thought to pivot on the existence of lo- 
cal icosahedral structures centered around solute atoms. 
Various idealized SRO configurations were presented in 
[46 . Schenk et al. experimentally verified icosahedral short 
range order (ISRO) in undercooled liquids [9 . Kelton et 
al. were the first to experimentally establish a connection 
between ISRO and the nucleation barrier [IQ . Later work 
further established the importance of ISRO in glasses [47l 

mm- 
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Many structural characterizations are oriented toward 
static viewpoint of the system, but some dynamical fea- 
tures have also been examined. Analysis of "free volume" 
(unoccupied space between atoms) fluctuations [45] have 
been used. Shear stress calculations investigate dynamical 
processes in glass forming materials [2T|[5Q] . Dynamic het- 
erogeneities involving cooperative motion of structures in 
a glass have also been studied [51^ 52. 53^ l54 [ l55j. 

Viable characterizations of SRO and MRO structures 
were advanced for low [56^ and high [49^ solute concentra- 
tions, binary [56,49 , and multicomponent systems [44 . 

Some methods of characterizing local structures in- 
clude Voronoi tesselation ^57|[3T| l56] , Honeycutt- Andersen 
indices [58^, and bond orientation ordering parameters 
[59] . These local measures center on an atom or a given 
link and, by definition, are restricted from detecting more 
complex longer range general structures. 

Experimental means to directly measure MRO struc- 
tures are given in [60^61 . Some potential MRO clusters 
were examined [22 ,56,49 . Some approaches to understand 
MRO use pattern matching to idealized MRO structures 
often constructed as agglomerations of perfectly ordered 
SRO features. 

A very useful experimentally driven approach for as- 
certaining MRO has been to look for "prepeaks" in X-ray 
and neutron scattering data [that is, peaks in the struc- 
ture factor S{q) that appear for wavenumbers corre- 
sponding to the inverse interatomic distances, which are 
lower in magnitude than that corresponding to the domi- 
nant S{q) peak]. We remark here that while this approach 
may capture general MRO structures, it is possible to have 
such structures without structure factor prepeaks (see Ap- 
pendix [C]). 



2.2 Preliminaries concerning our method 

Our unbiased structure characterization method extends 
multiresolution ideas [62] that have generally been re- 
served for network science applications (analyzing graphs 
in myriad social and biological networks) to complex ma- 
terials. Any complex physical system may be expressed 
as a network composed of nodes that code basic units 
of interest (e.g., atoms, electrons, etc.). Weighted links 
may capture the strength of the interactions between the 
different nodes or experimentally determined correlations 
(e.g., covariance or partial correlation contributions to the 
structure factor). After casting the system as a network, 
we then search for "communities" of nodes (e.g., clusters 
of atoms) that are more tightly linked to each other (or 
- in the case of the use of correlation functions as link 
weights - are more strongly correlated with one another) 
than to nodes in other clusters [63] . 

Our multiresolution method employs the notion of "com- 
munity detection", see , e.g., [64,65,66,67,68,63[70l|7K 
[721I731I741I75] . to quantitatively identify the "best" scale 
(or scales) for a complex physical system. Our approach 
does not rely on intuition or a knowledge of expected 
"important" features. Rather, it quantitatively estimates 



the best scale (s) through information-theory-based corre- 
lations, such as the variation of information (VI) [76^ or 
normalized mutual information (NMI), among different 
solutions. [In Appendix |Aj we review these information 
theory measures.] In essence, different copies of the com- 
munity detection problem are given to different solvers 
("replicas"). If many of these solvers strongly agree re- 
garding certain features of the solution, then these aspects 
are more likely to be correct. Similarly, a large discrep- 
ancy may indicate large fluctuations on a particular scale 
(or scales). Extrema in NMI or VI among the results of 
these independent solvers then indicate the best scales for 
the network. 

Multiple extrema can therefore indicate the existence 
of multiple relevant length/time scales. Although in most 
physical instances there is only a single dominant corre- 
lation length, there are many other cases in which more 
than one length/time scale is present [77] . 

One distinction between our work and some other es- 
tablished studies of local structures in glasses is that our 
analysis is not looking strictly at the positional structure. 
Rather, it evaluates structures in terms of the potential 
energies (i.e., the internal binding energies of the clusters, 
see also [22]). Our method can encapsulate weights that 
represent general statistical (pair of higher order) stress 
(or other) correlation functions, relative atomic displace- 
ments, etc. 

Our approach provides a perspective different from 
the "point-to-set" [78,79, 34 and other methods [41 . The 
point-to-set method examines the overlap between config- 
urations in a given volume (a "cavity" ) in an equilibrated 
system and compares those to configurations in the same 
cavity of the equilibrated system in which the boundary of 
the cavity was held fixed. Physically, it probes how prob- 
able it is to have a particular configuration within a disk 
or ball of a particular diameter given the boundary condi- 
tions. If many small clusters exist inside a sphere of fixed 
radius, then, a change in the boundary conditions will not 
significantly alter the bulk cluster distribution within the 
sphere. Conversely, if the sphere radius is smaller than 
the natural correlation length, then the number of config- 
urations compatible with the boundary will be small and 
the overlap will be large. A different but perhaps related 
approach, that of Ref. [41 , examines the distribution of 
structures inside a given volume to identify the correlation 
length. The method examines whether the distribution of 
configurations inside the volume occurs with a random 
frequency (when the linear scale of the volume is larger 
than the correlation length) or not (when the linear scale 
of the volume is smaller than the correlation length). 

Our method does not search for overlap at different 
scales for a multitude of configurations nor does it exam- 
ine their frequency. Rather, the pertinent structures are 
revealed by the information theory extrema between dif- 
ferent copies of the entire system. We do not tabulate 
possible configurations and their occurrence probabilities 
or examine the system in restricted volumes. 

Furthermore, the basic structures that we find may be 
used as the natural units in a renormalization group type 
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Fig. 1. A depiction of our simulated model glass former with 
three components "A", "B", and "C" with mixture ratios of 
88%, 7%, and 5%, respectively. The N = 1600 atoms are sim- 
ulated via IMD [82 in cube of approximately 31 A in size with 
periodic boundary conditions. The identities of the atoms are 
C (red), A (silver), B (green) in order of increasing diameters. 

analysis or coarse grained theories where clusters are re- 
placed by single nodes and an effective energy can be writ- 
ten that entails interactions between the different clusters 
alone. Finding the basic units is not trivial in amorphous 
systems such as glasses and numerous disordered systems. 
In these systems, there is generally no symmetry or any 
other obvious natural key that may determine how to op- 
timally partition the system on different scales. 
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Table 1. Fit parameters for Eq. (|lj) obtained from fitting con- 
figuration forces and energies to ab-initio data |83l[84j . The 
units of the parameters are such that given r in A, (/)(r) is in 
eV. (That is, the parameters ai, and as are dimensionless, 
ao is in A, a2 is in eVA^^ and as is in A~^.) The same-species 
(*) data is replaced by a suggested potential derived from gen- 
eralized pseudo-potential theory [85] (see also Appendix [P]) . 



tial energy function is [83] 

(/)(r) = [y) + ^ COS {asr + a^) , (1) 

where r is the distance between the centers of two atoms. 
This potential form incorporates a realistic weak long range 
interaction. Table [T] summarizes the parameter values 
which depend on the specific types for a pair of interacting 
atoms, and Fig. [2] shows the respective potential plots. 

The interaction parameters {a^}f^^ were determined 
[83] by fitting configuration forces and energies to ab- 
initio data [84]. The same-species model interactions are 
finally replaced by that suggested by generalized pseudo- 
potential theory (GPT) [85 . As illustrated in Fig. [l] we 
simulate N = 1600 atoms in a cubic system approximately 
31 A in size using periodic boundary conditions. This 
width is approximately twice the size of any suspected 
MRO structures. 

The system is initialized at a temperature of T = 1500 
K and allowed to equilibrate for a long time using a con- 



3 Simulations of model glasses 

We examine a model glass former derived from a three- 
component AlYFe metallic glass [80 which we designate 
as "A", "B", and "C" in mixture ratios of 88%, 7%, and 
5%, respectively. The presence of the different components 
B and C assists in the formation of a glassy state [81] since 
few pure compounds manifest a glassy state except under 
extreme preparation conditions. 



3.1 Ternary model glass former 

One system that we examine is derived from a three- 
component AlYFe metallic glass. As depicted in Fig. [l} 
it is a model glass former with components designated as 
"A", "B", and "C" in mixture ratios of 88%, 7%, and 5%, 
respectively. 

We use classical molecular dynamics (MD) [82 to sim- 
ulate the system dynamics. For this, we need accurate 
effective pair potentials that portray the pairwise interac- 
tions between the atoms in the system. Our model poten- 
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Fig. 2. The pair potentials for our three-component model 
glass former (see Fig.[T]). We indicate the atomic types by "A", 
"B", and "C" which are included with mixture ratios of 88%, 
7%, and 5%, respectively. The units are given for a specific 
candidate atomic realization (AlYFe) discussed in the text. 
The same-species data uses a suggested potential derived from 
generalized pseudo-potential theory [85] (see also AppendixlD]) . 
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stant number of atoms (N), a constant volume (V), and 
a constant energy (E). That is, we work within the NVE 
ensemble. After allowing for system equilibration, we save 
s high temperature configurations separated by a fixed pe- 
riod of simulation time. Prior to cooling, the length scales 
in the system are changed by 1% to account for the in- 
crease in density as a result of cooling since we choose 
to cool the system in an NVT ensemble to control the 
temperature. The system is then rapidly quenched to a 
temperature of T = 300 K, and it is allowed to equilibrate 
(in a mostly frozen state) in an NVE ensemble. We again 
save s separate low temperature configurations separated 
by a long period of simulation time. 



3.2 Lennard-Jones glass 

We additionally test the ubiquitous Lennard-Jones poten- 
tial using the Kob- Andersen (KA) 80:20 binary liquid [86 
which lies in the glass- forming mixture region [87 . The 
potential is 



12 



(2) 



where a oi j3 designate one of two atomic types A and B. 
Specifically, in accord with KA we set the dimensionless 
units CAA = 1.0, CAB = 0.50, cbb = 1-5, aAA = 1-0, 
aAB = 0.88, and asB = 0.80. 

As in the ternary glassy system above, we use MD [82j 
to simulate a LJ system of = 2000 atoms. The system 
is initialized at a temperature of T = 5 (using energy 
units where the Boltzmann constant /c^ = 1) and allowed 
to evolve for a long time. We save s high temperature 
configurations separated by 1000 time steps. The time step 
size is At = 0.0069 in LJ time units. Then, the system is 
rapidly quenched to a temperature of T = 0.01 which is 
well below the glass transition temperature of the KA-LJ 
system. The system is evolved in this mostly frozen state, 
and we save s low temperature configurations separated 
by 1000 steps of simulation time. 



4 Multiresolution clustering on amorphous 
materials 

Our idea is to apply, for the first time, multiresolution 
network analysis methods to ascertain all pertinent struc- 
tures in complex amorphous materials. A key subcompo- 
nent of this analysis is the community detection method 
itself. We first explain these ideas in network analysis and 
their physical analogs. 



4.1 Physical Motivation 

In an ideal decomposition of a large system into decoupled 
subsystems (communities), there is no interaction between 
different communities, and the system is effectively that 
of an ideal gas of disjoint communities. Stated differently. 



in the simplest setting in which the Hamiltonian would 
be block diagonal, the evolution of nodes (e.g., atoms) in 
each community would be decoupled from all other nodes 
in other communities. We next consider a fundamental 
Newtonian many particle setting. If the total force on a 
cluster is zero (i.e., if the cluster is strongly decoupled 
from all others) Fduster = then, the particles in that 
cluster will drift (on average) with a uniform velocity (the 
center of mass velocity of that cluster), i.e., v ^ const. 
[88 

In such instances, we may treat each community as a 
different particle in an ideal gas of non interacting such 
particles. The general problem is to find (the time de- 
pendent) permutation that renders the pair interaction 
strengths and/or correlations between particles into the 
best possible block diagonal form (on the time scale cho- 
sen). Community detection emulates this for graphs. In 
the atomic realization that we discuss in this article, this 
emulates a partition of the system into optimally decou- 
pled clusters such that the system may be viewed as an 
''ideal gas" of decoupled communities. 

For slow cooling of a liquid which enables crystalliza- 
tion, a first order or critical transition appears in the com- 
munity detection problem. A similar transition is materi- 
alized in many slowly cooled liquids as they crystallize to 
form nicely ordered solids. By contrast, in an extremely 
rapid cooling of a liquid, the interactions between the par- 
ticles are similar to those in spin-glass systems (i.e., the 
particles are not nicely organized and consequently the 
inter-particle interactions harbor a large degree of ran- 
domness). It is notable that a spin-glass transition appears 
in the community detection problem for random graphs 



4.2 Community detection 

In a graph such as that depicted in Fig. [3) nodes cor- 
respond to abstracted fundamental elements of the sys- 
tem, and edges represent defined relationships between 
the nodes. Community detection describes the problem of 
finding strongly connected groups of nodes. Nodes in dif- 
ferent clusters are more weakly connected than the nodes 
belonging to the same cluster. As alluded to earlier, we 
apply direct physical analogies where a node corresponds 
to a single atom. Edges, and their corresponding weights, 
are directly defined by the associated pair- wise interaction 
energy (or, in the absence of known interactions, may be 
set by measured inter-particle correlations). Specifically, 
we use the interatomic potential energies in Eqs. ([T]) and 
([2|. This potential model of the network edges is physi- 
cally appealing in that finding the best partition for the 
network is akin to minimizing the cluster binding energies 
of the physical system. 

Our Potts model Hamiltonian for community detection 
is [63] 



^({^}) = -\Y. i^^i - lBij)Siai,aj). 



(3) 
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Fig. 4. Panel (a) depicts our replica construction for the physical system at a "static" time to (no time separation between 
replicas, see Appendix [kJ . For presentation purposes, only a few nodes (atoms) are illustrated. Panel (b) depicts a similar set 
of replicas separated by a time At between successive replicas. In both panels (a) and (b), we generate a model network for 
each replica using the potential energy between the atoms as the respective edge weights for the network. Independently, within 
each replica, we then subsequently minimize Eq. ^ at a given value of 7 using the algorithm briefly described in Sect 



4.2 



Afterwards, we use the information measures in Appendix [A] to evaluate how strongly the set of replicas agree on the best 
partition. 




Fig. 3. An arbitrary weighted network with 4 natural commu- 
nities (strongly connected) depicted as distinct node shapes. 
The goal in community detection is to identify any such 
strongly related clusters of nodes based on their defined edge 
relationships. Solid lines depict weighted links corresponding 
to complimentary or attractive relationships where Aij > 
and Bij = in Eq. ([3|. Gray dashed lines depict missing, ad- 
versarial, or repulsive relationships where Aij = and Bij > 
m Eq. (|3|. In both cases, the relative link weight is indicated 
by the respective line thicknesses. For presentation purposes, 
missing intercommunity edges are not depicted. In this paper, 
we directly relate the edges (attractive and missing/repulsive) 
to the interaction energy between pairs of atoms which implies 
that the natural groups would correspond to bound clusters of 
atoms. 



Globally minimizing this Hamiltonian corresponds to iden- 
tifying strongly connected clusters of nodes. The elements 



of the matrices Aij and Bij are the edge weights and are 
defined as follows: an "attractive" weight has Aij > 1 if 

if the nodes are not 
> if the nodes 



nodes i and j are connected, Aij 
connected, a "repulsive" weight has Bij 
are not connected, and Bij = if nodes i and j are con- 
nected. In this paper, we define model the weights Aij and 
Bij by pair-wise potential energies (see Sect.l5|. In princi- 

icTud 



pie, we can generalize the Hamiltonian to incTude n-body 
correlations or interactions. 

We split the "attractive" (ferromagnetic) and "repul- 
sive" (anti-ferromagnetic) contributions into two separate 
weighted matrices in order to insert the model weight 7 
that adjusts the energy trade-off between the two types of 
interactions. The parameter 7 allows us to vary the tar- 
get scale of the community solution. The spin states 
designate the community membership of each node i with 
a range 1 < (Ji < q where q is the number of communi- 
ties. This number q may be variable [62""63' (such as in the 
multi-resolution that we perform in the current work) in 
order to find the optimal solutions or held fixed. Node i 
is a member of community k if ai = k. In this Hamilto- 
nian, each spin interacts only with other spins in its 
own community. 

Briefly, using Eq. (§ our community detection algo- 
rithm rapidly moves nodes between communities based on 
the current lowest energy assignment until no more moves 
are possible (see also Appendix IbI) . We then attempt t 
independent trials and select the lowest energy trial as 
the best division. [62 ,63 [The number of trials serves as 
our optimization parameter in this greedy algorithm. It 
effectively allows the algorithm to explore more possible 
configurations before selecting the best configuration for 
a given replica. The number of required trials in order 
to achieve a prescribed accuracy monitors the computa- 
tional complexity (correlating with the number of local 
minima in which individual trials may get stuck). Some- 
what better optimization could be obtained with a heat 
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bath solution algorithm [89 at a cost of substantially in- 
creased computational effort.] This community detection 
algorithm partitions the network into communities by as- 
signing a unique cluster membership for each node. Fur- 
ther details are provided in the appendices. 

Local features in metallic glasses generally exhibit in- 
terconnecting short range structures [56 . In our commu- 
nity detection problem, this feature corresponds to allow- 
ing "overlapping" node memberships where atoms can be 
members of more than one local cluster. We incorporate 
this effect by assigning a node as a secondary member of 
every community for which it has a negative binding en- 
ergy in terms our Potts model in Eq. (Is]) (see Appendix 

We can express the general partition function for a 
community partition with inter-community interactions 
which may correspond to the surface terms of clusters in 
Random First Order Transition theory (RFOT). Our pa- 
rameter 7 effectively plays the role of scaling the relation 
between surface and bulk terms in RFOT. A high value 
of 7 corresponds to large surface effects while a small 7 
corresponds to dominant bulk effects. 

Several strengths of our method are listed below. The 
analysis is independent of the type of structures that are 
being analyzed (structured, amorphous solid, and possi- 
bly even liquid systems) . It is robust to noise in the model 
network, and it yields very accurate results with a simple 
and fast greedy algorithm [63 . Because edge assignments 
are based on relative node positions (through the interac- 
tion potential), our method should be robust with respect 
to translational or rotational motion of solid structures in 
the system (such as crystal nucleation). 

4.3 Multi resolution community detection 

Multiresolution methods [9Qll91[l92ll62l[93j extend the ideas 
of community detection to identify the "best" division(s) 
over a range of network scales ("resolutions"). We test s 
independent realizations of the system ( "replicas" ) over all 
relevant network scales by specifying different values of 7 
in Eq. (|3|. Networks with a weakly defined structure will 
result in more diverse solutions among the different repli- 
cas. Conversely, strongly defined structures will result in 
replica solutions that agree more strongly. When replicas 
represent time-separated configurations, a strong agree- 
ment among the replicas corresponds to consistent phys- 
ical structures over time. The time separation between 
different snapshot can be tuned to find the corresponding 
pertinent structures for general time scales. The particu- 
lar case of vanishing time separation corresponds to dif- 
ferent representations of the same static system. See Fig. 
|4j In the particular case of vanishing time separation, we 
permute, for identical nodes (atoms), the numbers label- 
ing the nodes within the network; this leads to replicas 
corresponding to different i nitia l starting points for the 
algorithm described in Sect. |4.2[ [Another possibility now 
shown in Fig. |4]is that of replicas being different copies of 
the same system as it evolves in space-time (i.e., a replica 
not corresponding to a given snapshot of the system at 
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Fig. 6. Panels (a) and (b) show the plots of information mea- 
sures In, V, if, and / and the number of clusters q (right-offset 
axes) versus the Potts model weight 7 in Eq. The ternary 
model system contains 1600 atoms in a mixture of 88% type A, 
7% of type B, and 5% of type C with a simulation temperature 
oi T — 300 K which is well below the glass transition for this 
system. This system shows a strongly correlated set of replica 
partitions as evidenced by the information extrema at (i) in 
both panels. A set of sample clusters for the best resolution at 
7 ^ 0.001 is depicted in Fig.fTol 



a fixed time but rather monitoring the entire system as 
it evolves in space-time) . This latter possibility will be 
alluded to in Appendix [R[] 

Within this framework, we identify the best resolu- 
tion (s) by analyzing how well the independent replica so- 
lutions agree with each other in terms of information con- 
tent [62 . See Fig. [s] for a schematic involving time sepa- 
rated replicas. Extremal correlated resolution(s) identify 
the best division(s) of the network. We apply the vari- 
ation of information (VI) V metric and the normalized 
mutual information (NMI) In measure to evaluate the 
level of similarity among all pairs of the s replicas (see 
Appendix |A|. Extrema correspond to locally stable so- 
lutions that remain locally unchanged for variations in 
the system scale. In the case of very stable system resolu- 
tions, local extrema can be replaced by plateaus in these 
information theory measures that indicate no change in 
the system solution over an extended range of resolution 
scales 7 (as seen in the networks examined in [62 and, 
in the appendix of the current work, in some crisp lattice 
systems analyzed in Appendix [0|. We can further extract 
additional qualitative information about the "stability" of 
network partitions across a range of resolutions by exam- 
ining the average number of clusters q [ 9Q|[9T| l94] . mutual 
information / [90 , or the Shannon entropy H fQ^^^^. 
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Fig. 5. A depiction of our multiresolution algorithm using the replicas of Fig. [4|^b) for a range of resolution parameters 7 in 
Eq. ([3|. We solve the set of replicas at each 7^ after which we utilize the information measures in Appendix [A] such as In in 
the schematic, to measure the level of agreement among the replicas for each tested resolution. The NMI In or VI V extrema 
(or plateaus in some instances) indicate preferred (or more "stable") resolutions. 
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Fig. 7. Panels (a) and (b) show the plots of information mea- 
sures In, V, H, and / and the number of clusters q (right-offset 
axes) versus the Potts model weight 7 in Eq. ([3|. The ternary 
model system contains 1600 atoms in a mixture of 88% type A, 
7% of type B, and 5% of type C with a simulation temperature 
of T = 1500 K which is well above the glass transition for this 
system. At this temperature, there is no resolution where the 
replicas are strongly correlated. See Fig. [6] for the correspond- 
ing low temperature case where the replicas are much more 
highly correlated at 7 ^ 0.001. 



4.3.1 Phase diagram of the community detection problem 

Apart from trying to find the best possible division of com- 
munities in a particular system, we have also looked at the 
ease at which such a solution can be found. This further 
indicates when, physically, the system might be analyzed 




Fig. 8. A depiction of the full partitioned system where 
unique cluster memberships are depicted as distinct colors 
(best viewed in color). The atomic identities are B, A, C in or- 
der of increasing diameters. Overlapping nodes (multiple mem- 
berships per node) are added to these communities to deter- 
mine the best interlocking system clusters. 



in terms of nearly decoupled communities and when it 
cannot. It also enables us, by varying the parameters and 
temperature to map out when a solution is present or not, 
to determine whether the solutions found can be trusted 
and are physically relevant. 

When solving the system of Eq. ([3| at non-zero tem- 
perature for a given network (i.e., for an atomic configura- 
tion that is held fixed), entropic effects can, on their own, 
lead to a transition as the temperature is increased [89^. 
This transition corresponds to a spin glass type transi- 
tion for random systems capturing a transition from easy 
to hard computational problems. These transitions may 
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(b) 




Fig. 9. Panel (a) is the full system cube, and panels (b) - 
(d) show three sample clusters (one distinct cluster each using 
periodic boundary conditions) within the simulation box . Note 
that the algorithm can identify structures beyond immediate 
short range neighbors. 



transitions associated with the computational complexity 
of the community detection problem (as further manifest 
via the distribution of partitions g{A)) A decompo- 
sition of the type of Eq. Q is exact for a (standard uni- 
form) q state Potts model (JY = — ^ Si^j Aij^(Ti,aj where 
Aij = 0, 1) wherein the interaction energy between spins 
in different clusters (with each cluster /community given 
by uniform value of a^) is identically zero {Scn^aj = 0)- 



In Fig. 11 , we plot a specific measure [62 of the compu- 
tational complexity as a function of the density of energet- 
ically attractive inter-community links Pout and tempera- 
ture T for a particular random graph. Within the solvable 
flat region (starting at low Pout and T) the system can be 
effectively decomposed into decoupled elements, i.e., the 
partition function satisfies Eq. Q with a well-defined set 
of partitions {yl}. This region is separated by ridges of 
high complexity in which the community detection prob- 
lem becomes exceedingly hard from the unsolvable region 
in which no sensible community detection occurs (weak 
thermal fluctuations aid the system in avoiding metastable 
states while large thermal fluctuations are detrimental). 
While the specific phase diagram boundaries were found 
with the Hamiltonian of our method (discussed next), the 
phase diagram changes little when other methods are used. 
In several simple cases the phase boundaries coincide with 
those of the known cases (e.g., the phase transition of the 
Ising model on the square lattice when investigated as a 
network with links representing the strength of the spin 
exchange) . 



also be standard equilibrium transitions for sufficiently 
regular systems - critical for several regular graphs (e.g., 
the square lattice viewed as a graph analyzed via a Potts 
model for < 4 communities) or first order {q > 4 in 
the example above); for weak disorder, the latter transi- 
tions may be rounded off and further display signatures of 
Griffiths type behavior [101\ We briefly elaborate on this 
idea below. As illustrated in the examples analyzed in this 
work, spatially increasing low temperature structures in a 
supercooled liquid as well as correlation lengths in sys- 
tems such as the Ising spin systems can be ascertained 
via our method. In a general decomposition of an inter- 
acting system into optimally decoupled groups of particles 
("communities"), the partition function is approximated 
as 

{A} c=l 

Here, Zc is the partition function as computed with the 
Hamiltonian of the entire system for the particles in com- 
munity c, g{A) is the frequency of obtaining a particular 
set of communities yl in a decomposition of the entire sys- 
tem into optimally disjoint communities, and qA is the to- 
tal number of communities in the partition A. Thus, phys- 
ical transitions or crossovers may, in this approximation, 
be related to transitions (or divergent scales) or crossovers 
in the communities found in A themselves and/or phase 




Fig. 10. A depiction of some of the best clusters of the low 
temperature (T = 300K) ternary system at the peak replica 
correlation at feature (i) in Fig. [g] These clusters include over- 
lapping node membership assignments where each node is re- 
quired to have an overall negative binding energy to the other 
nodes in the cluster. The atomic identities are C (red), A (sil- 
ver), B (green) in order of increasing diameters. 
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5 Multiresolution application to model glass 
formers 

We assign edges between the nodes (atoms) with the re- 
spective weights based on the empirical pair-potentials 
given by Eqs. ([T]) and ([2|. Specifically, we calculate the 
potential energy between each pair of nodes i and j 
in the system and then shift each value by a constant 
^0 to obtain = + (assuming that ^ as 
the distance between particles i and j tends to infinity 
(r oo)). The shift > is necessary for the commu- 
nity detection algorithm to properly partition the network 
of atoms since it provides an objective definition of which 
interatomic spacings are preferable for a well-defined clus- 
ter and which are preferred to be excluded from a cluster. 

In our particular application here, we calculate the av- 
erage potential energy of the system and set 0o = — ^avg- 
For use in Eq. (§, we define an edge with a weight Aij = 
—(t^ij between nodes i and j if < 0, and we weight 
any missing links (or "repulsive edges") by Bij = if 

> 0. We then solve both model systems over a large 
range of 7 using s = 12 replicas and t = 10 optimization 
trials per replica. 

While (j)o = — 0avg is an intuitive shift that accom- 
plishes the goal of an objective cluster definition here, it 
is not an appropriate shift for some problems. For exam- 
ple, using 00 = — 0avg turns out to be problematic in some 
cases for lattice models. In a general setting, we examine 
a continuum of potential shifts 0o and monitor extrema 
in the information theory measures as a function of both 
7 in Eq. ([3| and (f)o. 

In addition to the systems tested below, we applied 
the algorithm to various test cases including square, tri- 
angular, and cubic lattice structures. The algorithm is 
able to correctly identify the natural leading order scales 
(plaquettes and composites of plaquettes as "cascades" in 
the information theory correlations). Further testing in- 




volved two-dimensional defects (dislocations, interstitials, 
etc.) and domain walls in a lattice. Defects in triangular 
lattices occurred most frequently near cluster boundaries. 

We also tested static configurations for the ternary 
model glass system where each replica is a model of the 
same configuration. There we detected structures in both 
low and high temperatures where the high temperature 
"structures" are more fragile (that is, harder to solve in 
the clustering problem). This corresponds to identifying 
relevant transient features in a dense liquid. 



5.1 Ternary model glass results 

In Figs, [g] and [7| panels (a) and (b) show the information 
theory based correlations (averaged over all replica pairs 
as in [62 ) over a range of network resolutions. The lower 
temperature system at T = 300 K in Fig. [6ja) shows a 
peak NMI at (za) with a corresponding VI minimum at 
{ih). Fig. fsl depicts an example of the full system parti- 
tion. Fig. [9] shows some sample clusters within the sim- 
ulation bounding box at resolution parameter value of 
jbest — 0.001 where we include overlapping node mem- 
berships (the replicas correlations are calculated on parti- 
tions), and Fig. 10 depicts additional samples of the best 



clusters. The corresponding high temperature (T = 1500 
K) solutions have a much lower NMI at jbest — 0.001 indi- 
cating significantly worse agreement among replicas. That 
is, one would expect that the high temperature system 
T = 1500 K is in a liquid state, so any observed features 
are not dynamically stable across all replicas (snapshots of 
the system over time). At T = 300 K, the best structures 
have consistent cluster sizes that are exclusively MRO. 

The plateau regions for 7 > 10 are similar to the LJ 
plot in Fig. [T2j but in this system the NMI plateau is 
lower. In the high temperature case in Fig.[7| there are ad- 
ditional "almost-plateaus" for the range 0.001 ^ 7 ^ 0.1. 
These plateaus represent a region of structural transition, 
but we are not concerned with them because the replica 
correlations are very low. 



5.2 Binary Lennard-Jones glass results 



In Figs. 12 and 14, panels (a) and (b) show the data for the 
replica information correlations over a range of network 
resolutions. The lower temperature system at a tempera- 
ture of T = 5 (in units of kB = 1) in Fig. |6ja) shows a 
plateau in NMI at (ia) with a corresponding VI plateau at 
(ih) which are the local extrema (y = is a trivial solu- 



tion with only one cluster in this problem). Fig. 13 depicts 



Fig. 11. A measure of the computational complexity x of the 
community detection problem as a function of the temperature 
T and the density of inter- community links Pout- 



a sample of the best clusters, including overlapping node 
memberships, at resolution {i) for jbest — 10^. The cor- 
responding higher temperature solutions at jbest — 10^ 
(see Figs. 14 anc 15) have a lower NMI (indicating a weak 
agreement among replicas). The dependence number of 
replicas (See Fig. 15) required to achieve high accuracy 
underscores the faint agreement between contending solu- 
tions and the high temperature complexity of the prob- 
lem. Our identified structures for this LJ model system 
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are consistent in terms of the cluster sizes and are almost 
exclusively SRO configurations with simple adjunct-type 
atoms extending into the low end of MRO size structures. 



6 Conclusions 

Our method is a new and very general approach to deter- 
mine the natural multi-scale structures of complex phys- 
ical systems. We do not bias the expected configurations 
in any way. The required input is that of inter-particle 
interactions (or measured correlations as further detailed 
in Appendix |N|. Information theory extrema (including 
plateaus) between contending solutions give the different 
pertinent structures on all important length scales (lat- 
tice scales and correlation lengths) of the system in an 
unbiased unified way. 

Apart from benchmarks on general networks and on 
lattices and spin s ystems with and without defects (see 
Appendices 0|P|Q ), to illustrate the feasibility of this ap- 
proach, we focused in this work on structural glasses. The 
detection of structure in structural glasses is a heavily in- 
vestigated hard problem. By the use of our method, we 
identified consistent SRO or MRO structures at temper- 
atures below the glass transition in two different model 
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Fig. 12. Panels (a) and (b) show the plots of information 
measures In, V, and / and the number of clusters q (right- 
offset axes) versus the Potts model weig ht 7 in Eq. The 
LJ system contains 2000 atoms in a mixture of 80% type A 
and 20% type B (Kob- Andersen binary LJ system [86]) with 
a simulation temperature of T = 0.01 (energy units) which is 
well below the glass transition of Tc — 0.5 for this system. This 
system shows a somewhat strongly correlated set of replica 
partitions as evidenced by the information extrema at (m,b) 
in panels (a) and (b). A set of sample clusters for the best 
resolution at 7 = 10^ is depicted in Fig. 13 




Fig. 13. Several of the best clusters for the peak replica cor- 
relation at feature (i) in Fig. 
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These clusters include over- 
lapping node membership assignments where each node is re- 
quired to have a overall negative binding energy to the other 
nodes in the cluster. The atomic identities are B (silver) and 
A (red) in order of increasing diameters. 



glass formers. Our analysis evaluates structures in terms 
the potential energies (i.e., the internal binding energies 
of the clusters). This approach differs from some other 
methods of structural analysis that look strictly at the 
relative atomic positions. We briefly comment on the re- 
lation between our work and that of simpler mode analy- 
sis of the interactions - the latter have been indeed been 
used to detect a link between spatial structure and lo- 
cal dynamics [53| [95]. Local heterogeneous dynamics was 
shown to be correlated with topological defects [96] and 
thus (as a consequence of [53^95 ) with the system modes. 
(Indeed, a recent work [97 directly reaflirms such a con- 
nection.) Similarly, some old variants of graph partition- 
ing methods such as direct spectral clustering [98""99" as 
well as community detection employ a normal mode anal- 
ysis [71]. Indeed, a method based on particle dynamics in 
high dimensions enabled multi-scale community detection 
[72] , Other different yet conceptually related approaches 
also include oscillator synchronization analogies [100,73 . 
In the above approaches, physical analogies were made. It 
is thus natural to suspect that the approach may inverted 
and that community detection will link spatial structure 
with dynamics in a broad class of physical systems where 
the physics based analysis of the interacting multi-particle 
system is hampered by the shear complexity of the prob- 
lem. The current work indeed fleshes out this link in de- 
tail and introduces a method for the direct detection of 
general spatio-temporal structures in rather general phys- 
ical systems. Notably, via the use of information theory 
correlations and extrema as a function of the Hamilto- 
nian parameters therein (7 of our deflning Hamiltonian 
of Eq. [3] and 0o of Sect, [s]), we are able to identify in an 
unbiased way all of the natural scales of the system. In 
more rudimentary approaches this needs to introduced by 
hand (e.g., a cutoff on mode occupancy in spectral based 
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approaches) and it is not obvious how to determine all 
pertinent structures of a complex physical system. 

Our approach identifies MRO as the dominant feature 
of our ternary model glass former with no strongly de- 
fined SRO. In contrast, the LJ system shows a largely SRO 
structure with adjunct atoms that create near-MRO struc- 
tures. Compounding the changes in structure that we find 
by analyzing the atomic system at different temperatures 
and minimizing the energy function to determine the op- 
timal division into clusters, there are also entropic effects. 
The distribution of optimal partitions becomes wider and 
less pronounced due to these effects as the temperature 
increases. This is a validation of our arguments in Sect. 

imi 

On a lattice, plateaus in information theory correlation 
steps correspond to a cascade of structures starting from 
the smallest dyads of nodes, to basic plaquette structures 
(square, triangle, etc.), and growing ever larger (two joined 
plaquettes etc.). In Ising spin systems at different temper- 
atures on a square lattice, the domains of and " 
spins are separated from one another by domain walls. The 
information theory plateaus correspond similarly to the 
cascade of small plaquette structures found on the lattice 
itself (i.e., the single plaquette, two joined plaquettes etc.) 
up to a cutoff scale set by the domain wall. This is sen- 
sible since no clear structure is found beyond the domain 
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Fig. 14. Panels (a) and (b) show the plots of information 
measures In, V, and / and the number of clusters q (right- 
offset axes) versus the Potts model weig ht 7 in Eq. The 
LJ system contains 2000 atoms in a mixture of 80% type A 
and 20% type B (Kob- Andersen binary LJ system 86 ) with a 
simulation temperature of T = 5 (energy units) which is well 
above the glass transition of Tc — 0.5 for this system. At this 
temperature, the replicas are significantly less correlated than 
the corresponding low temperature case in Fig. |12| 
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Fig. 15. A sample depiction of dispersed clusters for the LJ 
system Eq. ([2| at a temperature of T = 5 (in units where 
ks = 1)- T he s hown clusters correspond to the mult iresolut ion 
plot in Fig. [m] at value of the resolution parameter of 7 = 10^. 
These clusters are a sample of high temperature counterparts 
to the low temperature clusters in Figs. 12 and |13| Panels (a) 
and (b) show a more typical example of dispersed clusters at 
a number of replicas s = 10. In some cases, the identified high 
temperature clusters can be more compact, but not densely 
packed. Panels (c) and (d) provide sample solutions for s = 20 
replicas. An increasing replica (s) number (and generally also 
trial (t) number, see text) required to achieve better solutions 
is indicative of a greater computational complexity of the sys- 
tem. Physically, a larger time is required for the system to 
realize better clusters under ideal conditions (in the absence 
of quenching and any relaxation time constraints). The spar- 
sity of the identified clusters in this high temperature system is 
generally consistent across all clusters in the network solution. 
(This lies in contrast to the more compact and more strongly 
correlated structures found at low temperatures.) The atomic 
identities are B (silver) and A (red) in order of increasing di- 
ameters. 



length scale. The largest fluctuations occur at the bound- 
aries between different domains. These domain walls are 
directly attained by the extrema (those corresponding to 
the maximum in VI). Physically, they correspond to the 
scales at which the largest fluctuations occur where the 
large fluctuations lead to poor information theory correla- 
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Fig. 16. Plot of NMI In as a function of the potential shift 
00 and the Potts model weight 7 for the ternary model system 
in Sect. |5.l| The temperature is T = 300 K which is below the 
glass transition temperature for this system. This plot shows 
that the peak is roughly constant across a range of potential 
shifts. In general systems, we can obtain all natural scales by 
shifting both 00 and 7 and looking for extrema (and plateaus 
in some cases, see Sect.[0]). See Fig. [g] for the corresponding 

2D plot using 00 = (l)avg' 



tions between the different replicas. Fig. [42] corresponds to 
a sample depiction of the system at the maximum VI. The 
correlation lengths - set by the scale of the largest domain 
walls - are thus related to the scale which the fluctuations 
in the information theory overlaps are maximal (as man- 
ifest via maxima in VI). 

As we detail in Appendix [S| a general rigid amorphous 
solid supports shear and a divergent "shear penetration 
length". This length scale monotonically increases as a 
liquid is cooled to form a rigid glass. Within our graph 
theoretic method, we may employ the shear stress cor- 
relations as to represent graph weights to analyze such 
behavior. 
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A Information theory measures 

We utilize the variation of information [76 and normal- 
ized mutual information to measure the strength of the 



correlations among the independent replicas in our mul- 
tiresolution algorithm. The mutual information /(A, B) 
between community divisions A and B is 



QA QB 

/(AS) = ^^!^ log 

i=i j=i 



N 



(5) 



qa and qb are the number of communities in divisions A 
and 5, Tii and rij are the number of nodes in communities 
i and j, and TT/^j IS a "confusion matrix" that identifies 
the number of nodes in community i of partition A that 
are found in community j of partition B. For a single 
community division A, we can determine the self mutual 
information by H{A) = /(A, A) which is identically equal 
to the Shannon entropy. We use base 2 logarithms. 

We calculate the variation of information V{A^ B) by 



V{A, B) = H{A) + H{B) - 2/(A, B). 



(6) 



The range for VI is < V{A, B) < log A^. The normalized 
mutual information In{A^B) is 



In{A,B) 



2I{A,B) 
H{A)^H{B)' 



(7) 



The range for NMI is < In {A, B) <1. The NMI and VI 
afford slightly different perspectives on the replica corre- 
lations in the mult iresolut ion analysis above. 

These information measures are based on the cluster 
definitions and are not directly related to a thermody- 
namic entropy. The replica correlation measures could be 
improved by incorporating information about the assigned 
"overlapping" configurations (nodes may belong to more 




Fig. 17. Plot of NMI In as a function of the time (MD time 
steps) between configurations ts and the Potts model weight 7 
for the ternary model system in Sect. |5.l| The temperature is 
T = 1500 K which is above the glass transition temperature for 
this system. Intuitively, this plot shows that the correlations 
become weaker as the time between configurations is increased. 
This process of examining the correlations as a function of ts 
and 7 allows us to examine the relevant time scales in addition 
to the natural spatial scales that we identify in Sect. [5] 
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Fig. 18. Plot of NMI In as a function of the time (units of 10 
MD time steps) between configurations ts and the Potts model 
weight 7 for the ternary model system in Sect. |5.1| The tem- 
perature is T — 0.01 which is above the glass transition tem- 
perature for this system. Intuitively, this plot shows that the 
correlations become weaker as the time between configurations 
is increased. Intuitively, this plot shows that the correlations 
become weaker as the time between configurations is increased. 



than one cluster) such as in Ref. [92], by modifying the 
measures to account for the fact that same-type atoms 
are indistinguishable in this type of a model, or by utiliz- 
ing a thermodynamic entropy since we are dealing with a 
physical system in the current application. 



s 




Fig. 19. Plot of NMI In as a function of the number of trials 
s used to optimize the solution and the Potts model weight 7 
for the binary LJ system in Sect. |5.1| The temperatures are 
T — 0.01 in panel (a) and T = 5 in panel (b). The number 
of trials has a higher effect on the accuracy in the T = 5 
system, but the number of trials does not result in a drastic 
improvement to the accuracy of our algorithm for either LJ 
system. 



B Physical community detection models 



Minimizing the Hamiltonian of Eq. ([3]) using the algorithm 
briefly explained in Sect. |4.2[ models a dynamic commu- 
nity detection process where we search for a local energy 
minimum that indicates a "good" community partition 
(in general). In particular, we have applied a potential 
energy (PE) model of network edges, and we could, in 
principle, apply other edge definitions to obtain other rel- 
evant configurations with the following caveat. When we 
minimize the Hamiltonian using the edge weights, a low 
energy state corresponds to a good partition, so for con- 
sistency with the community detection problem, any edge 
definitions should ideally be extremized at the most fa- 
vorable configuration (s). The subtle distinctions between 
a few different natural edge weight models can result in 
different calculated clusters (beyond natural fluctuation 
caused by a high configurational entropy). 

An inverted PE model fits the above criterion, and it 
would at first appear to be conceptually consistent with 
a model based on the relative pairwise squared node dis- 
placement (SND) model. For both cases, small deflections 
about the minimum still indicate a good, even if not op- 
timal, configuration corresponding to the intuitive notion 
of a bound cluster. However, the SND model would cause 
a perfect crystal to be identified as a single contiguous 
cluster. The PE would not have this effect since more dis- 
tant nodes have a much lower potential energy. In effect, 
the PE model could identify the smaller scale units in a 
perfect crystal that the SND model could not isolate. 

Another intuitive edge model for bound clusters is re- 
lated to the attractive force that one may expect to ex- 
ist between constituent elements of a physical cluster. Of 
course, the system forces still cause the formation of a 
crystalline ground under the appropriate conditions, but 
the forces dictate the system's physical dynamics over 
time where the PE model (more) directly indicates the rel- 
evant ground statefs). Forces with a minimum such as the 
LJ model in Eq. ^ would be zero at the optimal ground 
state configuration (ideal crystal state), and the maximum 
force would be at a larger radius r than the PE minimum. 
While a force model of edge weights is also intuitive, it 
would correspond to (perhaps slightly) different clusters 
in a physical system compared to clusters derived from a 
PE model because the community detection Hamiltonian 
is extremized at different radii. 

This discussion of the subtle differences between the 
edge models leads to the natural question of which is the 
"best" physical cluster model. Another perspective is that 
the different edge weight models simply answer different 
questions. We have selected the PE definition since it best 
corresponds to the ideal ground state of the system in 
our pairwise interaction model. We then infer the best 
local clusters within configurations in local (frustrated) 
equilibrium (in the solid systems). 
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C MRO and structure factor pre-peaks 

One experimental approach to ascertain MRO is to look 
for "prepeaks" in the scattering data. That is, one can 
look for lower amplitude peaks in the structure factor S{q) 
which precede the dominant S{q) peak for wavenumbers q. 
While the approach may capture general MRO structures, 
it is possible to have MRO structures without significant 
pre-peak(s) in the structure factor plot. 

We illustrate the basic premise of this statement with 
an elementary example - that of a random arrangement 
of crystallites. Even though the Fourier transform of the 
mass density in each individual crystallite has sharp peaks 
at the reciprocal lattice vectors corresponding to these 
small crystallites, the structure factor obtained from the 
entire system may have those peaks vanish. We may, for 
instance, denote the locations of the centers of mass of 
individual grains i by Ri and denote the location of indi- 
vidual atoms j in each grain with respect to its center of 
mass by rj. In that case, the structure factor is 

Within each individual grain on its own, 

rjei 

will be sharply peaked about the corresponding reciprocal 
lattice vectors of grain i. However, the complete sum, 

S{k) =^eyip{ik-Ri)Si{k), (10) 

Ri 

may vanish if the relative distances between the different 
crystallites are random and lead to phases ex.p{ik • Ri) 
which obliterate the signatures of order in the individual 

S^{k). 



D Overlapping nodes between different 
communities 

We wish to account for the possibility of a given atom 
being connected to more than one physical cluster. For 
example, in a cubic lattice, each atom participates in the 
local structure of multiple unit cells. In community detec- 
tion, this corresponds to allowing "overlapping" commu- 
nity memberships where a node can be a member of more 
than one community. To accomplish this task, we select 
the lowest energy partition at the best resolution(s) of the 
model network [i.e., value(s) of 7 in Eq. (|3| corresponding 
to extrema in In or V]. 

First, we fix the initial node memberships including the 
number of communities q. We then sequentially iterate 
through the community memberships of each node and 
make changes according to the following: (1) place the 
node in any additional (non-member) clusters to which it 



is bound (a negative energy contribution), or (2) remove 
the node from any member clusters (except for the original 
membership) in which the current net energy contribution 
is positive. This process iterates through all nodes as many 
times as necessary until no node additions or removals are 
found. The total computational cost is slightly higher than 
the initial partitioning cost in Sect. [42] due to the multiple 
allowed memberships. See also [92j for another method 
that allows overlapping multiscale network analysis in a 
general fashion. 



E Potential shift in the ternary metallic glass 



In Fig. 16, we show the NMI correlations for a range of 
potential shifts ^0 and the Potts model weights 7 using 
s = 12 replicas and t = 10 trials for all data. The NMI 
peaks are roughly constant over a range of (j)o for the 
ternary model glass up to ^0 = 0-1 (the A- A interaction 
minimum 0min — 0.24 eV). 



F Time correlations in the ternary model 
glass 

In Fig. [T7| we show the NMI correlations for a range of 
separation times ts between replicas (units are in MD time 
steps) and the Potts model weights 7 using s = 12 replicas 
and t = 10 trials for all data. This plot intuitively shows 
that the correlations weaken as we increase the separation 




Fig. 20. Some of the best clusters for the peak V for the low 



temperature system shown in Fig. 12 in panel (b). The peak 
in the variation of information V generally correlates with the 
scale on which the fluctuations in the system division are most 
prominent. These clusters include overlapping node member- 
ship assignments where each node is required to have a overall 
negative binding energy to the other nodes in the cluster. The 
atomic identities are B (silver) and A (red) in order of increas- 
ing diameters. The resulting configurations constitute tightly 
bound objects. Representative corresponding high temperature 
clusters are shown in Fig. |21| 
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H An optimization effect in the binary LJ 
glass 



Fig. 21. The best clusters found for the peak V for the replica 
correlation within the high temperature LJ system (see text). 
These clusters include overlapping node membership assign- 
ments where each node is required to have a overall negative 
binding energy to the other nodes in the cluster. The peak in 
the variation of information V generally correlates with the 
scale on which fluctuations are the largest. The atomic identi- 
ties are B (silver) and A (red) in order of increasing diameters. 
The diffuse objects found at high temperatures are no longer 
as compact as at lower temperatures. The normalized mutual 
information at the peak V is also correspondingly lower than 
that for the low temperature system. 



time between replicas. Through varying the time step be- 
tween replicas, we examine the correlations as a function 
of ts and 7 which allows us to determine the most impor- 
tant time scale(s) (or relative equivalence). See Fig. |6] for 
the corresponding 2D plot using ts = 1000, and the static 
limiting case is shown explicitly in Sect. [K] 



G Time correlations in the binary LJ glass 

In Fig. [Tsj we show the NMI correlations for a range of 
separation times tg between replicas (units are in MD time 
steps) and the Potts model weights 7 using s = 12 repli- 
cas and t = 10 trials for all data. This plot shows that 
the correlations weaken only a little as we increase the 
separation time between replicas up to 10 000 MD time 
steps. Through varying the time step between replicas, we 
examine the correlations as a function of tg and 7 which 
allows us to determine the most important time scale (s) 
(or relative equivalence). See Fig. [g] for the corresponding 
2D plot using t = 1000, and the static limiting case is 
shown explicitly in Sect. [Kj 



In Fig. 19, we show the NMI correlations for a range of 
optimization trials t for s = 12 replicas and the Potts 
model weights 7. The plots show that the number of trials 
t has a small effect on the overall accuracy of the solution 
for either temperature (T = 0.01 or 5), but the effect is 
slightly more pronounced in the higher temperature T = 5 
case. See Figs. [12] and [13] for the corresponding 2D plots 
using s = 10. 

The high temperature clusters are generally much more 
dispersed. See Fig. 15 for sample clusters using t = 10 and 
20 at T = 5 where the clusters correspond to the multires- 
olution plot in Fig. [m] at 7 = 10^. The corresponding low 
temperature clusters are analyzed Fi g. [12] and presented 
in Fig. 13 Panels (a) and (b) in Fig. 15 display the typi- 
cal case of dispersed clusters at t = 10. In some instances, 
the high temperature clusters can be more compact, al- 
beit not densely packed, where panels (c) and (d) show 
two examples at 5 = 20. The clusters that are identified 
are generally consistent in terms of sparseness across all 
clusters in the solution. The frequency of occurrence of the 
more compact high temperature clusters, and its depen- 
dence on the level of optimization, is a subject for further 
study. 




(a) 



(b) 
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Fig. 22. Plots of the variation of information V between 
replica pairs as a function of the resolution parameter 7 for 
several system sizes in a LJ simulation (see Sect. 3.2). Panels 
(a), (b), and (c) use N = 2000, N = 4000, and N = 8000, re- 
spectively. Note that the value of 7 corresponding to the peak 
V scales downward with the system size (larger structures) 
which may indicate that it is correlated to a diverging length 
scale. 
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Fig. 23. Panels (a) and (b) show the plots of information mea- 
sures In, V, H, and / and the number of clusters q (right-offset 
axes) versus the Potts model weight 7 in Eq. ([3|. The ternary 
model system contains 1600 atoms in a mixture of 88% type 
A, 7% of type B, and 5% of type C with a simulation temper- 
ature of T = 300 K which is well below the glass transition for 
this system. In this system, we use a single static "snapshot" of 
the system to analyze what our multiresolution algorithm finds 
at T = 300 K. This low temperature case shows a preferred 
resolution at low 7 as evidenced by the information extrema 



at (i) in both panels. See Fig. 24 for the corresponding high 
temperature case. 
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I Multiresolution VI maxima in the LJ system 

In Fig. [T2j we show the multiresolution correlations as a 
function of the Potts model weight 7 for a LJ system of 
= 2000 nodes using s = 12 replicas and t = 10 trials 
for all data. The value of 7 corresponding to the peak in 
V corresponds to the "maximum complexity" of the en- 
ergy landscape which is often correlated to the system size 
(see Appendix Fig. 20 shows a sample of the best clus- 
ters corresponding to the VI peak which are roughly 300 
atoms in size (7 — 8 atoms in diameter) which is approxi- 
mately 1/7 of the size of the system. Fig. 21 shows some 
corresponding high temperature T = 5 clusters which we 



note are approximately twice the size (n 600 nodes) and 
more dispersed. 

As an additional note, we remark that, intuitively, one 
would expect that the VI maxima are directly related to 
the NMI minimum in this region of 7. This is not the 
case in this region because the VI metric does to zero as 
the network collapses to one cluster. It is rapidly shift- 
ing from a region of maximum complexity at a peak VI 
value, where most of the information comes from the sum 
of the Shannon entropies H{A) and H{B)^ for two par- 
titions A and 5, to minimum complexity at a value of 
V = 0. NMI has a different behavior here, as the ratio 
In = 2I{A,B)/[H{A)^H{B)], because the mutual in- 
formation between different replicas becomes smaller and 
approaches zero as the system collapses into larger com- 
munities. The low temperature LJ system in Fig. 12 shows 
a transitional very low peak in NMI which corresponds to 
a near bisection of the system. However, in this case the 
overlapping cluster configurations collapse to the entire 
system. 
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Fig. 24. Panels (a) and (b) show the plots of information mea- 
sures In, V, and / and the number of clusters q (right-offset 
axes) versus the Potts model weight 7 in Eq. The ternary 
model system contains 1600 atoms in a mixture of 88% type 
A, 7% of type B, and 5% of type C with a simulation tempera- 
ture of T = 1500 K which is well above the glass transition for 
this system. In this system, we use a single static "snapshot" 
of the system to analyze what our multiresolution algorithm 
finds at T = 1500 K. This high temperature case does show 
a preferred resolution at low 7 as evidenced by the informa- 
tion extrema at (i) in both panels, but the NMI maxima are 
somewhat lower. The solutions also take longer to solve ac- 
curately meaning that the community detection energy land- 
scape is more complicated. See Fig. [23] for the corresponding 
low temperature case. 
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Fig. 25. Panels (a) and (b) show the plots of information 
measures In, V, and / and the number of clusters q (right- 
offset axes) versus the Potts model weight 7 in Eq. ([3|. As in 
Sect. |5.2| this LJ system contains 2000 atoms in a mixture of 
80% type A and 20% type B with a simulation temperature of 
T = 0.01 (energy units) which is well below the glass transition 
for this system. In this test, we use a single static "snapshot" of 
the system. This low temperature system has a higher overall 
level of correlation among the replicas than the high tempera- 
ture case in Fig. |25| and the overall results are very similar to 
what was observed for the time separated replicas in Fig. [12] in 
Sect. 15.2] See Fig.|27|for a sample of some of the best observed 
clusters. 
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Fig. 26. Panels (a) and (b) show the plots of information 
measures In, V, and / and the number of clusters q (right- 
offset axes) versus the Potts model weight 7 in Eq. ([3|. The LJ 
system contains 2000 atoms in a mixture of 80% type A and 
20% type B with a simulation temperature of T = 5 (energy 
units) which is well above the glass transition for this system. 
In this analysis, we use a single static "snapshot" of the system. 
This high temperature system has a lower overall correlation 
than the low temperature case in Fig. |25| and the overall results 
are very similar to what was observed for the time separated 
replicas in Fig. [l4]of Sect. |5.2| 



J Finite size effects in the LJ system 



In Fig.|22j we show the VI correlations as a function of the 
Potts model weight 7 for several system sizes of = 2000, 
N = 4000, and N = 8000 in panels (a), (b), and (c), 
respectively. We also used 5 = 12 replicas and t = 10 
trials for all data. The value of 7 corresponding to the 
peak in V scales downward with the system size (larger 
structures) indicating that the maximum fluctuations are 
scaling with the size of the simulation box. In the limit of 
an infinite size simulation box, the VI peak would tend to 
7 = 0. 



K Static Multi-resolution analysis on ternary 
model glass 

The inter-replica NMI/VI correlations change with the 
time separation interval between the replicas where the 
longer time separations intuitively result in poorer correla- 
tions. The other limiting case is for a static analysis of the 
system. Thus, we are interested in determining what our 
multiresolution analysis finds when we examine a single- 
time snapshot of the system. A single snapshot may, po- 
tentially, capture a transient feature of the system. This 



is illustrated in panel (a) of Fig. [4] in which the time sepa- 
ration between different replicas is set to zero. We applied 
the same algorithm as in Sect. [4] except that all replicas 
correspond to the same system time, and we solved the 
systems with s = S replicas and t = 4 trials per replica. In 
Figs. [23l and [24} our analysis identifies a static length asso- 
ciated with the instantaneous configuration. At both high 
and low temperatures, the static multi-resolution analysis 
display information extrema at low 7 (large spatial scales). 
The higher temperature system displays weaker correla- 
tions with sparser structures whereas the structures of the 
low temperature system display more significant correla- 
tions (and more compact structures). As seen in Figs. 23 
and [24] the high temperature system displays weaker cor- 
relations in the vicinity of the peak NMI than those of the 
low temperature system. That is, the high temperature 
solutions require more effort to solve accurately (there are 
far more metastable minima and more trials are required 
in order to find better contending solutions). (A similar 
occurrence was found for the high temperature LJ system 
in Fig. [TsJ) The solved clusters (not depicted) are compa- 
rable to those identified for the time-separated replicas in 
Sect. [5] for both system temperatures. 
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Fig. 27. A depiction of some of the best clusters for the peak 
replica correlation in the static system of Fig. |25| These clusters 
include overlapping node membership assignments where each 
node is required to have an overall negative binding energy 
to the other nodes in the cluster. The atomic identities are A 
(red) and B (silver). 



tentials is shown in Fig. 30 These interactions conform to 
those 

The lower temperature system at T = 300 K in Fig. 
shows a peak NMI at {ia) 
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aj snows a peak IN Ml at (zaj w ith a corresponding VI 
minimum at {ih). Following Sect. |5.1[ Fig. 31 depicts a 
sample of the best clusters at jbest — 0.001 where we in- 
clude overlapping node memberships (the replicas correla- 
tions are calculated on partitions). We used 5 = 12 replicas 
and t = 20 trials per replica. The corresponding T = 1500 
K high temperature solutions have a much lower NMI at 
jbest — 0.001 indicating very poor agreement among the 
replicas. At T = 300 K, the best structures have consis- 
tent cluster sizes that are MRO or a little larger which 
are generally larger than 1/2 the simulated system width. 
Therefore, it would beneficial to test the consistency of 
the clusters in a larger system requiring a substantially 
longer computational time. 

In this system the NMI plateau at 7 > 100 is actually 
higher than the configuration at (z), but the clusters are 
almost exclusively small (n 5 nodes) and are not com- 
pletely contiguous. The distinction in the results between 
the different potential models is likely due to the longer 
range Al-Al minimum in the ab initio potential parameter 
fit data as compared to the GPT |85j minimum. 



L Static Multi-resolution analysis on LJ glass 

As in Appendix [K] for the ternary system, we further 
tested a static version of the binary LJ model where the 



results are similar to the dynamic case in Sect. 5.2 both 
in the MRA plot and the resulting clusters. A notable ex- 
ception is that the distinction between the high and low 
temperature cases becomes more subtle in the multireso- 
lution plots. 

We applied the same algorithm as in Sect. [4] except 
that all replicas correspond to the same system time, and 
we solved the systems with s = 12 replicas and t = 10 tri- 
als per replica. In Figs. 25 and[26j we show our multires- 
olution analysis associated with the instantaneous config- 
urations. These two static plots are similar, but the high 
temperature case has slightly weaker correlations (lower 
NMI peak and higher VI minimum). In Figs. |25| and [26 



In Figs. 25 ^ 

where we use time separated replicas, the contrast in the 
MRA plots at different temperatures is stronger. Fig. [27 
shows a sample of the best clusters in the low tempera- 
ture case which are comparable to those identified for the 



time- separated replicas in Fig. 13 in Sect. 5.2 



M An analysis of yet another ternary system 



In F igs. [28] and [291 we repeat the analysis of Sects. [371] and 
5.1 , for an alternate ternary model glass system (AlYFe). 
In particular, we use the best parameter fit for the same- 
species interactions that are tabulated in Table |2] as op- 
posed to implementing the GPT [85 as used in Figs.[6]- 
Tq| and tabulated in Table [T] A plot of the alternate po- 
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Fig. 28. Panels (a) and (b) show the plots of information mea- 
sures /at, V, and / and the number of clusters q (right-offset 
axes) versus the Potts model weight 7 in Eq. ([3|. The ternary 
model system contains 1600 atoms in a mixture of 88% type 
A, 7% of type B, and 5% of type C with a simulation tem- 
perature of T = 300 K which is well below the glass transition 
for this system. This alternate system uses the best fit param- 
eter data for the same-species interactions as opposed to the 
pseudo-potential interaction used in Figs. [6l- |lQ| This system 
shows a locally preferred resolution as evidenced by the infor- 
mation extrema at (i) in both panels. A set of sample clusters 
for the best resolution at 7 0.0001 is depicted in Fig.|3l] The 
region for 7 > 100 actually has a higher correlation, but the 
clusters are very small (n :^ 5 nodes) and somewhat dispersed. 
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Fig. 29. Panels (a) and (b) show the plots of information mea- 
sures /at, V, H, and / and the number of clusters q (right-offset 
axes) versus the Potts model weight 7 in Eq. ([3|. The ternary 
model system contains 1600 atoms in a mixture of 88% type 
A, 7% of type B, and 5% of type C with a simulation temper- 
ature of T = 1500 K which is well above the glass transition 
for this system. This alternate system uses the best parameter 
fit for the same-species interactions as opposed to the pseudo- 
potential interaction used in Figs.[6l- |10| At this temperature, 
there is no resolution where the replicas are strongly correlated. 
See Fig. [28] for the corresponding low temperature case where 
the replicas are much more highly correlated at 7 :^ 0.001. 



N Multi-resolution analysis via measured pair 
correlation function 



ordering are of extreme interest. Atomic configurations 
that are consistent with the experimentally determined 
scattering data for liquid Zr8oPt2o at 250 K below the liq- 
uidus were generated using conventional Reverse Monte 
Carlo (RMC) methods [lOSi llOQllnQj . The result of anal- 
ysis uses unconstrained RMC which may or may not em- 
ulate the precise microscopic structure. We then analyzed 
a representative system with N = 7500 nodes using the 
algorithm outlined in Sect. |4.3[ 

Specifically, we analyze the system as a static model. 
That is, each replica is based on the one particular system 
representation since the data is obtained by RMC meth- 
ods rather than by a dynamical simulation. The multi- 
resolution analysis (MRA) is seen in Fig. 32 A few of 
the relatively uniform large clusters are shown in Fig. 33 
where each cluster is approximately 700 — 800 atoms in 
size. Interestingly, the ZrPt system shows a well-defined 
secondary MRA NMI peak near 7 :^ 1. Although the peak 
NMI value has a poor overall correlation among the replica 
solutions (low /at), it is notable that the system displays 
this secondary peak since it is entirely absent in the LJ 
binary liquid results in Sect. |5.2| 



O Multiresolution application to lattice 
systems 

We define several uniform lattices systems for the pur- 
pose of comparing the results to the model glasses where 
we use relatively small systems for presentation purposes. 
We would normally just treat respective model networks 
as unweighted networks, but here we wish to maintain a 
consistent analysis across all systems in the paper, so we 
further apply a "potential" shift (j)Q which corresponds to 



As stated in Sect. |2j we may apply the same multireso- 
lution network clustering ideas for other structure mod- 
els. In particular, in Fig. [32| we apply the method to ex- 
perimentally adduced pair correlation functions gahi'f) in 
an amorphous ZrPt system between the different of the 
different components a and h (Le., Zr-Zr, Zr-Pt, and Pt- 
Pt). That is, we set in Eq. (|3|, Aij = —{gij + Qij) if 
{Qij + 9ij) < and Bij = {gij - g^j) if {gij - g^j) > 0. 
Here, gij denotes the pair correlation between the atoms 
corresponding to nodes i and j, gij is a background av- 
erage (which we set to zero). Some sample clusters then 
found by our method are seen in Fig. [33] Due to the large 
system size, we used 5 = 8 replicas and t = 4 trials per 
replica. 

Amorphous systems with extended atomic order be- 
yond the nearest neighbor shell provide an excellent frame- 
work to test the identification of natural structural ele- 
ments. The ZrPt is a system that has been shown to have 
MRO in both the glassy and liquid state [103 , 104 , lOSlfTOGl 
1107) . MRO, or correlations beyond direct chemical bond- 
ing manifest as pre-peaks in the static structure factor 
extracted from scattering studies of liquids and glasses. 
The dominant structural elements at the heart of such an 
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Fig. 30. A plot of the model potentials for our three- 
component model glass former (see Fig. [T]) using the fit data in 
Table [2] We indicate the atomic types by "A", "B", and "C" 
which are included with mixture ratios of 88%, 7%, and 5%, 
respectively. The units are given for a specific candidate atomic 
realization (AlYFe) discussed in the text. Here, we apply the 
ab initio fit data for the same-species interactions, as opposed 
to the GPT i85j model from Table [l] in Fig.|2] 
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Fig. 31. A depiction of some of the best clusters for the peak 
replica correlation at feature (i) in Fig. [g] These clusters in- 
clude overlapping node membership assignments where each 
node is required to have an overall negative binding energy 
to the other nodes in the cluster. The atomic identities are 
C (red), A (silver), B (green) in order of increasing diame- 
ters. These clusters are generally larger than 1/2 the simulated 
system width; therefore, it would beneficial to test their con- 
sistency in a larger simulation (requiring substantially longer 
computational time). 




Fig. 33. A sample of clusters determined by a multiresolution 
clustering analysis based on the pair correlation function in 
Fig. [32] The system is a Zr8oPt2o system at 250 K below the 
liquidus (see text) where Zr atoms are depicted as red, and Ft 
atoms are depicted as green. Each cluster has approximately 
700 - 800 atoms. 
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Fig. 32. A plot for the multiresolution analysis of a binary 
Zr8oPt2o system at 250 K below the liquidus. Atomic configu- 
rations were generated using conventional RMC methods that 
are consistent with the experimentally determined scattering 
data for liquid Zr8oPt2o at 250 K below the liquidus. The plot 
shows a poorly correlated, but nevertheless well-defined, peak 
in NMI near 7^1. 



the negative of the average weight over all pairs of nodes 
(any non- neighbor has a weight of Bij = 1). 

We allow zero energy moves for the lattice solutions 
and perform a more strenuous optimization using s = 20 
replicas and t = 40 trials. These representative networks 
result in "imperfect" tilings of the favored local structures 
due to the constraints imposed by the perfect symmetry 



in the Hamiltonian and by the local solution dynamics. 
In the depictions, different colors represent distinct clus- 
ters (best viewed in color) and edges between clusters are 
made partially transparent as a visualization aid. No over- 
lapping nodes are assigned in the lattice depictions. 



0.1 Square lattice 

We define a uniform, initially unweighted, square lattice 
with N = 400 nodes. Edges are assigned to each neighbor 
in the x and y directions with periodic boundary con- 
ditions. The "potential" shift is 0o = 0.97995. We then 
perform the same multiresolution analysis to the graph as 
in the previous systems. In Fig. |34j we see that where are 
three dominant plateaus in the information measures. 

Except for the plateaus, the overall multiresolution 
pattern resembles the L J system in Fig. |5.2[ except for the 
presence of the plateaus in the lattice plot, which suggests 
that the LJ system is "more ordered" than the metallic 
glass model. However, a purely random graph [62 also 
shows a similar pattern, including a plateau, which in- 
dicates that there may exist an analogy between purely 
"random" and perfectly "ordered" systems in our anal- 
ysis. However, these data are not alone sufficient to be 
conclusive. 

We select a configuration at 7 = 60 corresponding to 
the center plateau. The lattice is depicted in Fig. [35] with 
q = 120 clusters of 78 squares, 4 triads, and 38 dyads. 
At this resolution, the square dominates the configuration 
which shows that our algorithm is able to isolate the basic 
unit cells of the lattice in a natural fashion. 

The plateau for 7 > 100 corresponds to essentially all 
dyads, and the plateau for 7 ^ 30 corresponds to a mix- 
ture of dyads, squares, and tight six- node configurations 
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(a square plus two adjacent nodes). The lower 7 plateau 
favors the six-node configuration in terms of the cluster 
energies, but the larger features are even more difficult 
with which to tile the lattice than for squares. 
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Fig. 34. The multiresolution analysis of a square lattice with 
periodic boundary conditions. Neighbors have an initial weight 
of Aij = 1 and non-neighbors have an initial weight of Bij = 
1. We further apply a "potential" shift of 00 = 0.97995 in 
order to be consistent with our previous analysis on glasses. 
There are three distinct plateaus in the information measures. 
A depiction of the system at 7 = 60 for the center plateau is 
shown in Fig. |35| 




Fig. 35. A partition of a square lattice with periodic bound- 
ary conditions. The corresponding multiresolution plot is seen 
in Fig. [34] We use the algorithm described in the paper at 
7 = 60 to solve the system. To aid in visualization of the clus- 
ters, neighbor links not in the same cluster are made partially 
transparent. In this configuration, there were ^ = 120 clusters 
with 78 squares, 4 triads, and 38 dyads which indicates that 
square configuration dominates the partition, and it shows how 
our algorithm can naturally identify the basic unit cells of the 
square lattice. 



forming its own single-node cluster. Together, the different 
plateaus show that our algorithm is able to isolate the 
basic unit cells of the lattice in a natural fashion. 



0.2 Triangular lattice 

Similar to the square lattice we define a uniform, ini- 
tially unweighted, triangular lattice with N = 400 nodes. 
Edges are assigned to each triangular neighbor using pe- 
riodic boundary conditions. The "potential" shift is ^0 = 
0.969925. We again perform the same multiresolution anal- 



ysis to the graph with the results shown in Fig. 36 There 
are three dominant plateaus in the information measures; 
and the overall multiresolution pattern resembles both the 
square lattice and the LJ systems except for the presence 
of the plateaus here. 

We select a configuration at 7 = 50 corresponding to 
the "left" plateau. The lattice is depicted in Fig. [37| with 
q = 104 clusters of 17 triads, and 86 "diamonds", and 1 
five-node cluster. The five node cluster is a result of a pref- 
erence being given for an isolated node to join a diamond 
as opposed to forming its own single- node cluster. At this 
resolution, the diamond configuration dominates. 

Two plateaus to the right of 7 > 65 are both strongly 
dominated by triads of nodes. The distinction between 
the two is that the central plateau favors an isolated node 
joining a triangle, to form a rare diamond, rather than 



0.3 Cubic lattice 



We further define a uniform, initially unweighted, cubic 
lattice with N = 1000 nodes. Edges are assigned to each 
neighbor in the x, and z directions using periodic bound- 
ary conditions. The "potential" shift is ^0 = 0.987 988. We 
perform the same multiresolution analysis to the graph as 
in the previous systems where the results are summarized 
in Figs. |38] and 1^ 

Out of q = 177 clusters, we identified 66 squares, 51 
six-node configurations (square plus two adjacent nodes), 
45 cubes, and 15 other assorted configurations smaller 
than cubes. At 7 50, the cube is the preferred clus- 
ter in terms of the cluster energy, but they consist of 
only slightly more than 25% of the clusters because the 
large cube configuration is more difficult to identify due 
to the perfect network symmetry and local constraints im- 
posed by the evolution of the community structure dur- 
ing the algorithm dynamics. The "middle" plateau repre- 
sents a square-dominated region (201 out of 294 clusters 
are squares), and the "right" plateau consists of dyads of 
neighbor nodes almost exclusively. 
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Fig. 36. The multiresolution analysis of a triangular lattice 
with periodic boundary conditions. Neighbors have an initial 
weight of Aij — 1 and non-neighbors have an initial weight of 
Bij = 1. We further apply a "potential" shift of 00 = 0.969925 
in order to be consistent with our previous analysis on glasses. 
There are three distinct plateaus in the information measures, 
but the latter two are closely related (see text). A depiction of 
the system at 7 = 50 for the left plateau is shown in Fig. |37| 




Fig. 37. A partition of a triangular lattice with periodic 
boundary conditions. The corresponding multiresolution plot 
is seen in Fig. |36| We use the algorithm described in the pa- 
per at 7 = 50 (the "left" plateau) to solve the system. To aid 
in visualization of the clusters, neighbor links not in the same 
cluster are made partially transparent. In this configuration, 
there were q = 104 clusters with 17 triads, 86 "diamonds," 
and 1 five-node configuration (see text) which indicates that 
the diamond configuration dominates the partition. Our algo- 
rithm can naturally identify the different basic unit cells of 
the lattice. By varying 00 (and the corresponding 7) we detect 
natural structures on varying scales. 



P Multiresolution application to a 2D Ising 
lattice 

We define two-dimensional square lattice of Ising spins 
using the Hamiltonian 



H 



(11) 



for all pairs of neighbor spins i and j. Neighbors with the 
same spin have an initial weight of Aij = 1, neighbors with 
opposite spins have an initial weight of Bij = 1, and non- 
neighbors have an initial weight of Bij = 0. There is an 
unavoidable inherent discontinuity in defining the inter- 
actions for the Ising system (only neighbors can interact) 
and the continuous systems elsewhere in this paper {all 
pairs of nodes interact). With this caveat, we analyze this 
Ising system according to the same multiresolution anal- 
ysis, and we then compare the results to the other tested 
systems. We apply s = 8 replicas and t = 8 trials. 

We could effectively treat the above Ising interactions 
using an unweighted network model since we can factor 
our the weight, but we again wish to maintain a consis- 
tent analysis across all systems in this paper. With this 
in mind, we further apply the average "potential" shift 
(j)o which varies on each defined network and which corre- 
sponds to the negative of the average weight over all pairs 
of nodes (including non-neighbors). We allow zero energy 
moves for the solution dynamics. No overlapping nodes 
are assigned in the lattice. 



In Fig.|40j there are two main plateaus for 7 > 100 and 
a noticeable configuration "shift" at 7 100. The "right" 
plateau for 7 > 500 consists largely of same-spin dyads 
except where a given spin has no matching neighbors. The 
"middle" plateau for 100 ^ 7 ^ 500 corresponds to a 
"natural" grouping of medium size clusters. In Fig. [4l] 
different colors represent distinct clusters (best viewed in 
color) and edges between clusters (neighbor spins with the 
same sign but that are in different clusters) are depicted 
in gray. Missing edges are not shown 

Here, the plateaus correspond to a cascade of struc- 
tures starting from the smallest dyads of nodes, to basic 
plaquette structures (square, triangle, etc.), and growing 
ever larger (two joined plaquettes etc.). In Ising spin sys- 
tems at different temperatures on a square lattice, the do- 
mains of "+" and " spins are physically separated from 
one another by domain walls. The plateaus here similarly 
correspond to the cascade of small plaquettes found on 
the lattice up to a cutoff scale set by the domain walls. 

The domain walls are closely related to the maximum 
in VI (NMI displays a different, non-extremal, behavior 
where the standard deviation is large). Physically, they 
correspond to the scales at which the largest fluctuations 
occur where the large fluctuations result in poor infor- 



mation correlations between the replicas. Fig. 42 shows 
a sample depiction of the system at this VI peak. Corre- 
lation lengths are thus likely related to poorly correlated 
replicas best indicated by a VI maximum. 
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Fig. 38. The multiresolution analysis of a cubic lattice with 
periodic boundary conditions. Neighbors have an initial weight 
of Aij — 1 and non- neighbors have an initial weight of Bij — 1. 
We further apply a "potential" shift of 00 = 0.987988 in order 
to be consistent with our previous analysis on glasses. There 
are three distinct plateaus in the information measures. The 
leftmost short plateau is the cube preferred resolution (in terms 
of cluster energy) with 45 out of 177 clusters are cubic clusters 
(no configurations larger than cubes are found). A depiction of 
the system at 7 = 50 for the left plateau is shown in Fig. [39] 



Q Multiresolution analysis of Lennard-Jones 
systems with elastic defects 



We define a uniform lattice based on a two-dimensional 
Lennard-Jones system using the LJ potential in Eq. ([2| 
where the lattice is constructed as the ideal ground state 
of the system using periodic boundary conditions. The 
main purpose is to check for consistency of our multires- 
olution method, but the lattice also allows us to perform 
a preliminary investigation of the behavior of our analysis 
in how it places lattice defects that are randomly inserted 
into the system. We used 5 = 8 replicas and t = 8 trials. 

The ground state of the LJ system is given by a tri- 
angular lattice with a lattice constant that is equal to the 
distance that minimizes the LJ interaction between two 
atoms. We define a uniform triangular lattice with N = 
2943 nodes where we create random "defects" by remov- 
ing some nodes. Edges are again assigned and weighted 
according to the LJ potential in Eq. ([2| applying periodic 
boundary conditions. (The "potential" shift is, as before, 
set to (j)o = 0.969925.) We perform the multiresolution 
analysis on the graph where the results are shown in Fig. 

m 

Next, we select a configuration at 7 31.6 from the 
peak corresponding to the largest natural communities for 
this system at this specific lattice spacing. An example of 
the best lattice partition is depicted in Fig. l44j Defects ap- 




Fig. 39. A partition of a cubic lattice with periodic bound- 
ary conditions. The corresponding multiresolution plot is seen 
in Fig. [38] We use the algorithm described in the paper at 
7 = 50 (the short "left" plateau) to solve the system. To aid 
in visualization of the clusters, neighbor links not in the same 
cluster are made partially transparent. In this configuration, 
there were q = 177 clusters with 66 squares, 51 six-node config- 
urations, 45 cubes, and 15 other assorted configurations. The 
cubic cluster is the preferred partition (in terms of cluster en- 
ergy) , but it is difficult to fill the system with cubes in practice 
due to the perfect symmetry of the network. 



pear to occur most likely near the boundary of neighboring 
communities. 



R Multiresolution analysis in space-time 

In addition to the potential interaction energies, we may 
model our multiresolution community detection analysis 
using the mechanical action 



1 



(12) 



Here, Mi is the mass of particle i and ri{t) denotes its 
location as a function of time. As throughout, (/) is the 
two particle interaction. Below, we refer to the case of 
a system in D spatial dimensions. In a discretized form 
(allowing, in particular, for discrete time steps At along 
the time axis), the action 



ai3 



(13) 



a/3 



with the nodes a and /3 representing nodes in space-time 
and Wais an effective "interaction weight" between them. 
Thus, we represent our physical system as a network em- 
bedded in {D + 1) dimensional space-time. We extremize 
this action "S" corresponding to the Hamiltonian of Eq. 
([3|. The spatial links are determined by (j){ri^rj) simi- 
lar to before (apart from a sign inversion) while the links 
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Fig. 40. A plot for the multiresolution analysis of a square 
lattice of Ising spins with periodic boundary conditions at a 
simulation temperature of T = 2.269 (close to the transition 
temperature). Neighbors with the same spin have an initial 
weight of Aij = 1 , neighbors with opposite spins have an initial 
weight of Bij — 1, and non-neighbors have an initial weight of 
Bij = 0. There are two distinct plateaus in the information 
measures for 7 > 100 with a significant configuration shift 
near 7 ^ 100. A depiction of the system at 7 = 400 for the 
center plateau is shown in Fig. [41] 



parallel to the time axis are determined by the particle 
masses. As before, extrema of the information theory cor- 
relation pinpoint the natural structures and scales in the 
system. In this case, there is only one tuning parameter 
7 as applied to different replicas of complete space-time 
networks. Thus, the space and time scales cannot be in- 
dependently adjusted and are intertwined and are related 
by a ratio which is the typical speed of sound c (or an 
appropriate average of such speeds) in equilibrated solids. 



S Continuum elasticity about local minima 
and the shear penetration depth 

The focus of our article is on a detailed bona fide descrip- 
tion of amorphous materials that invokes graph theory 
methods to ascertain general structure. The various exam- 
ples shown have hopefully clearly outlined the strengths 
of our approach. In this section, we depart from this ap- 
proach and outline a continuum formalism. 

Continuum elasticity in regular solids does not assume 
detailed knowledge of the underlying atomic constituents. 
In this section, we follow suit and outline general consider- 
ations for amorphous solids. In an ideal crystalline system, 
any small deformation about the crystalline ground state 
will raise the energy in a harmonic manner. The increase 
of energy (for non-plastic deformations) is captured by 




Fig. 41. A depiction of a partition of a square lattice of Ising 
spins with periodic boundary conditions at a simulation tem- 
perature of T = 2.269 (proximate to the critical temperature). 
The corresponding multiresolution plot is seen in Fig. [40] We 
use the algorithm described in the paper at 7 = 400 (cor- 
responding to the transition to the plateau with the highest 
NMI / lowest VI value) to solve the system. For presenta- 
tion purposes, we depict all spins by spheres and links that 
cross a community boundary as grey edges. As with the multi- 
resolution analysis of the square lattice, this multiresolution 
analysis shows that the dominant communities are square pla- 
quettes within same-spin domains. 



the elastic constants of the material. Nothing prevents, in 
principle, the application of these considerations for de- 
formation about local energy minima - as these pertain 
to an amorphous solid in a local energy minimum (an in- 
herent structure as it would have pertained to supercooled 
liquids). In what briefly follows, we consider the harmonic 
expansion of Eq. (12) about such a local energy minimum 



state. We then review general considerations about "shear 
penetration depth"- the length scale on which the medium 
reacts to shear. Earlier treatments considered the shear 
penetration depth in homogeneous media [llipil2j . Here, 
we consider what may occur when replicating these con- 
siderations for deformations about a general non-uniform 
local ground state. 

In what follows, we label the displacement of the i-th 
atom, Ui about its position in a given state, Ri by 



Ri 



(14) 



In Eq. ( |14[ ), denotes the location of the displaced atoms. 
Henceforth, we consider a coarse grained description in 
which n is a field defined in the continuum. (That is, we 
replace the lattice indices i in Eq.(14) by continuous spa- 
tial coordinates x.) The displacement vector u is di D- 
dimensional vector in real space. We apply field theoretic 
ideas introduced by [113j. In particular, we next follow 
the procedure of [30 rrT|[ll2j and write the correspond- 
ing action of the system in a manner similar to that of the 
energy of an anisotropic solid in {D -\- 1) dimensions. We 
will label the time direction as the a = direction and 
denote all spatial Cartesian coordinates hy a = 1,2,...,!). 
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Fig. 42. A partition of a square lattice of Ising spins with 
periodic boundary conditions at a simulation temperature of 
T = 2.269. The corresponding multiresolution plot is seen 
in Fig. |40| We use the algorithm described in the paper at 
7 = 3.60 X lO"'^, corresponding to the peak VI, to solve the sys- 
tem. For presentation purposes, we depict all spins by spheres 
and links that cross a community boundary as grey edges. The 
identified clusters vary substantially across the replicas (max- 
imum fluctuations) , but they are of the scale of the size of the 
spin domains. This qualitatively implies that the poor informa- 
tion correlations are naturally correlated with the individual 
domains of and "-" spins. At this scale, the fluctuations 
are largest (and hence the related weakest information theory 
correlations) . This scale on which fluctuations are maximal (co- 
inciding with the Ising domain walls) should be contrasted with 
that of Fig. [41] for which the fluctuations were minimal (low- 
est variation of information (VI)) and the normalized mutual 
information (NMI) was the highest. 



We will use Latin indices a = 1,2,...,!) to label the spa- 
tial directions alone. In the continuum limit, expanding 
the energy to harmonic order about any local minima (in- 
herent state), the action reads [ 111( 130]. 



(15) 



In Eq. (15), repeated indices are to be summed over. The 



quantities Caf3ab derive from effective elastic moduli. In 
regular solids, C^^ah are space (and time) independent 
quantities. In an amorphous medium, however, Ca(3ab{x^ t) 
are generally space (and time) dependent. The effective 
elastic constant Cooa6 in Eq . ([l5|) is set by the mass den- 
sity, i.e., Cooab = p{x^t). 130(11 lH[ll2j When a and (3 
both assume values between 1 and D, the quantities Caf3ab 
are set by the elastic moduli of the solid (e.g., the shear 
modulus jii). The square root of the typical ratio between 
the usual spatial elastic moduli and those along the time 
axis for the space-time representation of the solid in Eq. 
( p!5| is given by a typical speed of sound in the system 

One approach to defining defects such as dislocation 
and disclinations in regular solids is to consider the net 
change in the displacement or local orientation associated 



with closed trajectories. Defect densities relative to a given 
inherent state may be defined in a similar way. More for- 
mally, we can dualize (via a Hubbard Stratonovich trans- 
formation) the quadratic action of Eq.(15) and associate 
the dual fields with defect current densities in space- time. 
Replicating the analysis performed for regular solids in 
130(1111( 1112 to amorphous solids with an invert ible space- 
time dependent elastic moduli CajSabi we obtain the same 
result obtained for regular solids. Namely, dislocation and 
disclination current densities are given by 



ja 



with the local rotation given by 



(16) 



(17) 



The temporal direction {a = 0) components of these ten- 
sors are the usual (dislocation and disclination) defect 
densities of elasticity. 



0? 



(18) 



The condition of a divergence-less stress tensor (a) in the 
absence of an applied force, enables us to write it as a curl 
of a gauge field. 



gf3;^Xa 



(19) 



Elastic shear is mediated by the gauge field B. When a 
plastic component of the displacement field u exists, it 
couples minimally to the stress-gauge fields via 



X dt B°"'Ja, 



(20) 



for the coupling between the dislocation currents and gauges. 
Mathematically, the structure is very similar to that in 
electromagnetism when currents couples minimally to gauge 
fields. Just as charges screen applied electromagnetic fields 
(as, e.g., in the Debye screening of plasmas), defect charges 
can screen the elastic shear. In the uniform medium, when 
a defect condensate (llZ^'ol 7^ 0) appears, it gives rise a 
screening of the elastic shear |3Q|llllj , with a screening 
length set by 



Ct 



(21) 



In Eq. (21 ), ct is the transverse sound velocity and /i is the 
shear modulus. In the case of simple dilation of the elas- 
tic moduli when Ca(3ab{x^t) = w{x^t)Kai3ab where Ka/s 
is a space-time independent tensor and is a dilation 
function, the calculation of the elastic shear penetration 
depth can be replicated for the amorphous system to yield 
Eq. (21) when the quantities refer to the local values of 
As/iear, Ct, V^o, and /i. A similar effect appears for general 
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space-time dependent elastic moduli |114j when consider- 
ing elastic deformations about a local energy minimum. 
As is seen from Eq. (21), when the defect density tends to 



zero, the elastic shear penetration depth diverges 



X 



shear 



OO 



(22) 



as the system becomes rigid throughout. In order to as- 
certain detailed shear response in a system, we may em- 
ploy the shear stress correlation functions as link weights 
within our community detection algorithm. This is topic 
of a future study. 
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Fig. 43. A multiresolution analysis of a 2D triangular LJ lat- 
tice with periodic boundary conditions. Edges are weighted 
according the LJ potential in Eq. ^ (with a potential shift 
00 = 0.969925). There are two preferred regions: a small peak 
on the left and a large plateau on the right, where the peak here 
corresponds to the largest possible "natural" clusters. The de- 
fects make only a small alteration to the multiresolution plot. 
A depiction of the system at 7 ^ 31.6 for the left peak is shown 
in Fig. [44] 




Fig. 44. A depiction of a partition of a 2D LJ triangular 
lattice with periodic boundary conditions. The corresponding 
multiresolution plot is seen in Fig. |43| We use the algorithm 
described in the paper at 7 ^ 31.6 (the left peak) to solve 
the system. Our algorithm places defects generally near the 
boundary of the communities in order to minimize the energy 
cost of the defects in the community assignments. 
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